/* File: Functions.c */ 
/*
####????

*/

#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include <Rdefines.h>

//define the number of loci to use
enum{ n_disp_loci = 10 }; 
enum{ n_neutral_loci = 10 };

/*--------------------------structures----------------------------------------*/
/*define the structure of individual
	x & y are location & sex -- 0 is male & 1 is female
	Da & Db are alleles for dispersal linked loci
	Na & Nb are alleles for neutral loci */
typedef struct {
	int x, y;
	short sex;
	short Da[n_disp_loci], Db[n_disp_loci], Na[n_neutral_loci], Nb[n_neutral_loci];
} ind;

/* structure of the neighbor list... simply holding 40 x,y values */
typedef struct {
	int n, x[8], y[8];
} nb_list;

/*--------------------------function prototypes-------------------------------*/
//void init_inds (ind * cur_inds, int spX, int spY, int n);

/*--------------------------function definition-------------------------------*/

/* define the function to initialize a neighbor list given 8 next neighors for zones 
0	1	2
3	4	5
6	7	8 */
void initialize_neighborlist (nb_list *nb, int spX, int spY){
	nb[0].n = 5;
	nb[0].x[0]=0; nb[0].y[0]=1;
	nb[0].x[1]=0; nb[0].y[1]=spY-1;
	nb[0].x[2]=1; nb[0].y[2]=1;
	nb[0].x[3]=1; nb[0].y[3]=0;
	nb[0].x[4]=1; nb[0].y[4]=spY-1;
	nb[1].n = 8;
	nb[1].x[0]=-1; nb[1].y[0]=spY-1;
	nb[1].x[1]=-1; nb[1].y[1]=0;
	nb[1].x[2]=-1; nb[1].y[2]=1;
	nb[1].x[3]=0; nb[1].y[3]=spY-1;
	nb[1].x[4]=0; nb[1].y[4]=1;
	nb[1].x[5]=1; nb[1].y[5]=spY-1;
	nb[1].x[6]=1; nb[1].y[6]=0;
	nb[1].x[7]=1; nb[1].y[7]=1;
	nb[2].n = 5;
	nb[2].x[0]=-1; nb[2].y[0]=spY-1;
	nb[2].x[1]=-1; nb[2].y[1]=0;
	nb[2].x[2]=-1; nb[2].y[2]=1;
	nb[2].x[3]=0; nb[2].y[3]=spY-1;
	nb[2].x[4]=0; nb[2].y[4]=1;

	nb[3].n = 5;
	nb[3].x[0]=0; nb[3].y[0]=-1;
	nb[3].x[1]=0; nb[3].y[1]=1;
	nb[3].x[2]=1; nb[3].y[2]=-1;
	nb[3].x[3]=1; nb[3].y[3]=0;
	nb[3].x[4]=1; nb[3].y[4]=1;
	nb[4].n = 8;
	nb[4].x[0]=-1; nb[4].y[0]=-1;
	nb[4].x[1]=-1; nb[4].y[1]=0;
	nb[4].x[2]=-1; nb[4].y[2]=1;
	nb[4].x[3]=0; nb[4].y[3]=-1;
	nb[4].x[4]=0; nb[4].y[4]=1;
	nb[4].x[5]=1; nb[4].y[5]=-1;
	nb[4].x[6]=1; nb[4].y[6]=0;
	nb[4].x[7]=1; nb[4].y[7]=1;
	nb[5].n = 5;
	nb[5].x[0]=-1; nb[5].y[0]=-1;
	nb[5].x[1]=-1; nb[5].y[1]=0;
	nb[5].x[2]=-1; nb[5].y[2]=1;
	nb[5].x[3]=0; nb[5].y[3]=-1;
	nb[5].x[4]=0; nb[5].y[4]=1;

	nb[6].n = 5;
	nb[6].x[0]=0; nb[6].y[0]=-1;
	nb[6].x[1]=0; nb[6].y[1]=1-spY;
	nb[6].x[2]=1; nb[6].y[2]=-1;
	nb[6].x[3]=1; nb[6].y[3]=0;
	nb[6].x[4]=1; nb[6].y[4]=1-spY;
	nb[7].n = 8;
	nb[7].x[0]=-1; nb[7].y[0]=-1;
	nb[7].x[1]=-1; nb[7].y[1]=0;
	nb[7].x[2]=-1; nb[7].y[2]=1-spY;
	nb[7].x[3]=0; nb[7].y[3]=-1;
	nb[7].x[4]=0; nb[7].y[4]=1-spY;
	nb[7].x[5]=1; nb[7].y[5]=-1;
	nb[7].x[6]=1; nb[7].y[6]=0;
	nb[7].x[7]=1; nb[7].y[7]=1-spY;
	nb[8].n = 5;
	nb[8].x[0]=-1; nb[8].y[0]=-1;
	nb[8].x[1]=-1; nb[8].y[1]=0;
	nb[8].x[2]=-1; nb[8].y[2]=1-spY;
	nb[8].x[3]=0; nb[8].y[3]=-1;
	nb[8].x[4]=0; nb[8].y[4]=1-spY;
}

/* define the neighborhood zones (list above) of the grid data... */
void assign_nb_values (int **nb, int spX, int spY){
	for (int xx=0; xx<spX; xx++){
		for (int yy=0; yy<spY; yy++){
			if (xx==0){
				if (yy==0) { nb[xx][yy] = 0;		
				} else if (yy==spY-1) { nb[xx][yy] = 6;
				} else { nb[xx][yy] = 3; }
			} else if (xx==spX-1) {
				if (yy==0) { nb[xx][yy] = 2;		
				} else if (yy==spY-1) { nb[xx][yy] = 8;
				} else { nb[xx][yy] = 5; }				
			} else {
				if (yy==0) { nb[xx][yy] = 1;		
				} else if (yy==spY-1) { nb[xx][yy] = 7;
				} else { nb[xx][yy] = 4; }			
			}
		}
	}
}

/* check the defined neighborhood zones */
void check_nb_values (int **nb, int spX, int spY){
	printf("values in the neighbor zones are...\n");
	for (int yy=0; yy<spY; yy++){
		for (int xx=0; xx<spX; xx++){
			printf("%i ",nb[xx][yy]);
		} printf("\n");
	} printf("\n");
}

/*Initialises a matrix with n individuals
	individuals have a location (spX, spY) and are 
	initilised with gene for p.disp
	based on Stable range edge functions */
void init_inds (ind *individuals, int spX, int spY, int n) {
	int ii, jj; //some iteration variables
	for (ii = 0; ii < n; ii++) {
		individuals[ii].x = (int)(unif_rand() * spX) % spX;	individuals[ii].y = (int)(unif_rand() * spY) % spY; // get random x & x locations
		individuals[ii].sex = Rf_rbinom(1,0.5); //define the sex
		for (jj = 0; jj < n_disp_loci; jj++) { //assign the initial allelles
			individuals[ii].Da[jj] = Rf_rbinom(1,0.5); individuals[ii].Db[jj] = Rf_rbinom(1,0.5); 
		} 
		for (jj = 0; jj < n_neutral_loci; jj++) { //assign the initial allelles
			individuals[ii].Na[jj] = Rf_rbinom(1,0.5); individuals[ii].Nb[jj] = Rf_rbinom(1,0.5);
		}
	}
}

/* copy the new individuals over the current */
void copy_inds (ind *ci, ind *ni, int n) {
	for (int ii = 0; ii < n; ii++) {
		ci[ii].x = ni[ii].x; ci[ii].y = ni[ii].y ; ci[ii].sex = ni[ii].sex;
		for (int jj = 0; jj < n_disp_loci; jj++) { ci[ii].Da[jj] = ni[ii].Da[jj]; ci[ii].Db[jj] = ni[ii].Db[jj]; } 
		for (int jj = 0; jj < n_neutral_loci; jj++) { ci[ii].Na[jj] = ni[ii].Na[jj]; ci[ii].Nb[jj] = ni[ii].Nb[jj]; } 
	}
}

/* quick function to print out the individuals */
void check_inds (ind *individuals, int spX, int spY, int n) {
	int ii, jj; //some iteration variables
	for (ii = 0; ii < n; ii++) {
		Rprintf("individual %i -- X = %i -- Y = %i -- sex = %i\nlinked loci = ",ii,individuals[ii].x,individuals[ii].y,individuals[ii].sex);
		for (jj = 0; jj < n_disp_loci; jj++) { //assign the initial allelles
			Rprintf("%i%i",individuals[ii].Da[jj],individuals[ii].Db[jj]);
		} 
		Rprintf("\nneutral loci = ");
		for (jj = 0; jj < n_neutral_loci; jj++) { //assign the initial allelles
			Rprintf("%i%i",individuals[ii].Na[jj],individuals[ii].Nb[jj]);			
		}
		Rprintf("\n\n");		
	}
} 

/* reset the density to 0 */
void reset_grid(int **grid, int spX, int spY) {
	int xx, yy; //some iteration variables
	for (xx=0; xx<spX; xx++) {
		for (yy=0; yy<spY; yy++) { grid[xx][yy] = 0; }//populate with 0s 
	}
}

/* check the grid values */
void check_grid(int **grid, int spX, int spY) {
	printf("values in the grid are...\n");
	for (int yy=0; yy<spY; yy++){
		for (int xx=0; xx<spX; xx++){
			printf("%*i ",4,grid[xx][yy]); 
		} printf("\n");
	} printf("\n");
}

/* summarize the number of individuals and number of males */
void count_individuals(ind *individuals, int **dens, int **nm, int n) {
	for (int ii=0; ii<n; ii++) { //loop through all individuals & extract density info
		if (individuals[ii].sex>=0) { dens[individuals[ii].x][individuals[ii].y]++; } //add one to density if individual is alive
		if (individuals[ii].sex==0) { nm[individuals[ii].x][individuals[ii].y]++; } //add one to males if sex == 0
	}
}

/* function to return the first empty postion in mInds vector */
int first_empty(int *tt, int len){ 
	for (int ii=0; ii<len; ii++) {
		if (tt[ii] == -99) { return(ii); }
	}
	return(-99);
}

/* cycle through individuals and index positions of males */
void index_males(ind *individuals, int ***Index, int **nm, int n) {
	for (int ii=0; ii<n; ii++) { //loop through all individuals & index male positions
		if (individuals[ii].sex==0) { 
			int xx = individuals[ii].x; int yy = individuals[ii].y; //get the locations
			Index[xx][yy][first_empty(Index[xx][yy],nm[xx][yy])] = ii; 
		}
	}
}

/* check the index of the number of individuals */
void check_mIndex (ind *individuals, int ***Index, int **nm, int xx, int yy) {
	printf("sex for the individuals are...\nindividual    index    sex\n");
	for (int ii=0; ii<nm[xx][yy]; ii++) { printf("%i\t\t%i\t%i\n",ii,Index[xx][yy][ii],individuals[Index[xx][yy][ii]].sex); }
}

/* hassel commins population growth where b=1 (contest competition) generates expected fecundity */
double hass_comm(double lambda, int K, int N){ return lambda/(1+(lambda-1)*N/K); }

/* mutates a locus with probability m_rate*/
int mutate(int curr, double m_rate){
		if (curr==1) curr=0;
		if (curr==0) curr=1;
	return(curr);
}

/* function to define the characteristics of the new individuals given we now know how many offspring each female will have */
void reproduce(ind *ci, ind *ni, int *offspring, int n, int **nm, int ***Index, double m_rate) {
	int nn=0; //counter for ni
	int mutn=0;
	int d_mut_loc[2*n_disp_loci], n_mut_loc[2*n_neutral_loci];
	for (int ii=0; ii<n; ii++){
		if (offspring[ii]>0){ //females that have offspring only
			int male = (int)(unif_rand() * nm[ci[ii].x][ci[ii].y]) % nm[ci[ii].x][ci[ii].y]; //choose male
			male = Index[ci[ii].x][ci[ii].y][male];
			for (int jj=0; jj<offspring[ii]; jj++){ //step through offspring and allocate to nis
				ni[nn].x=ci[ii].x;
				ni[nn].y=ci[ii].y;
				ni[nn].sex=Rf_rbinom(1, 0.5);
				for (int ll=0; ll<n_disp_loci; ll++){ //assign dispersal alleles from both parents
					if (Rf_rbinom(1, 0.5)==1) ni[nn].Da[ll]=ci[male].Da[ll]; //Da=father's contribution (i.e., gamete)
					else ni[nn].Da[ll]=ci[male].Db[ll];
					if (Rf_rbinom(1, 0.5)==1) ni[nn].Db[ll]=ci[ii].Db[ll];	//Db=mother's contribution
					else ni[nn].Db[ll]=ci[ii].Da[ll];
				}
				for (int ll=0; ll<n_neutral_loci; ll++){ //assign neutral alleles from both parents
					if (Rf_rbinom(1, 0.5)==1) ni[nn].Na[ll]=ci[male].Na[ll]; //Na=father's contribution (i.e., gamete)
					else ni[nn].Na[ll]=ci[male].Nb[ll];
					if (Rf_rbinom(1, 0.5)==1) ni[nn].Nb[ll]=ci[ii].Nb[ll];	//Nb=mother's contribution
					else ni[nn].Nb[ll]=ci[ii].Na[ll];
				} 
				//Mutation
				mutn = Rf_rbinom(2*n_disp_loci, m_rate);
				if (mutn > 0){ //mutation at dispersal loci
					for (int ml=0; ml<(2*n_disp_loci); ml++){ //clear mutation location array
						d_mut_loc[ml] = -9;
					}
					for (int mm=0; mm<mutn; mm++){ //get locations for dispersal mutations
						d_mut_loc[mm] = (int)(unif_rand() * n_disp_loci) % n_disp_loci;
					}
					int cc=0;
					while (d_mut_loc[cc]>=0){ //mutate dispersal alleles
						if (Rf_rbinom(1, 0.5)==1) ni[nn].Da[d_mut_loc[cc]] = mutate(ni[nn].Da[d_mut_loc[cc]], m_rate);
						else ni[nn].Db[d_mut_loc[cc]] = mutate(ni[nn].Db[d_mut_loc[cc]], m_rate);
						cc++;
					}
				}
				mutn = Rf_rbinom(2*n_neutral_loci, m_rate);
				if (mutn > 0){  //mutation at neutral loci
					for (int ml=0; ml<(2* n_neutral_loci); ml++){ //clear mutation location array
						n_mut_loc[ml] = -9;
					}
					for (int mm=0; mm<mutn; mm++){ //get locations for neutral mutations
						n_mut_loc[mm] = (int)(unif_rand() * n_neutral_loci) % n_neutral_loci;
					}
					int cc=0;
					while (n_mut_loc[cc]>=0){ //mutate neutral alleles
						if (Rf_rbinom(1, 0.5)==1) ni[nn].Na[n_mut_loc[cc]] = mutate(ni[nn].Na[n_mut_loc[cc]], m_rate);
						else ni[nn].Nb[n_mut_loc[cc]] = mutate(ni[nn].Nb[n_mut_loc[cc]], m_rate);
						cc++;
					}
				}
				nn++; //increment ni counter
			}
		}
	}
}



/* disperse new individuals to neighboring cells */
void disperse(ind *ni, int n, nb_list *nb, int **nbI, double d_cost) {
	for (int ii=0; ii<n; ii++){ //loop through all individuals
		double dprob=0, dcost=0;
		for (int jj=0; jj<n_disp_loci; jj++){
			dprob+=ni[ii].Da[jj];
			dprob+=ni[ii].Db[jj];
		}
		dprob=0.3*(dprob/(2 * n_disp_loci)) + unif_rand() * 0.7; // dirty means of getting heritability of 0.3
		dprob=Rf_rbinom(1, dprob);
		dcost=Rf_rbinom(1, d_cost);
		if (dprob==1 && dcost==0){
			int tt = nbI[ni[ii].x][ni[ii].y];//get the type zone for selecting a neighbor
			int tt2 = (int)(unif_rand() * nb[tt].n) % nb[tt].n; // get a random neighbor
			ni[ii].x += nb[tt].x[tt2]; ni[ii].y += nb[tt].y[tt2]; // alter the current position
		}
		if (dprob==1 && dcost==1){ //marks dead individuals for removal later
			ni[ii].sex=-1;
		}
	}
}

/* function to create a matrix of key individual information */
/*
SEXP output_individuals(ind *ci, int n) {
	double *out; SEXP ans;
	int ncols = 2 + n_disp_loci + n_neutral_loci;
	PROTECT(ans = alloc3DArray(REALSXP, n, ncols, 2));
	out = REAL(ans); //pointer to output dataset
	for (int ii=0; ii<n; ii++) { //cycle through each of the individuals
		if (ci[ii].sex >=0){ //if the individual is alive, copy the necessary info
			for (int kk=0; kk<2;kk++){ //cycle through each chromosome
				out[ii + 0 * n + kk*ncols*n] = ci[ii].x; out[ii + 1 * n + kk*ncols*n] = ci[ii].y; // set the xy values
				if (kk==0){
					int tcol = 2; //define the column place holder
					for (int jj=0; jj<n_disp_loci; jj++) { out[ii + tcol * n + kk*ncols*n] = ci[ii].Da[jj]; tcol ++; } // get the loci information
					for (int jj=0; jj<n_neutral_loci; jj++) { out[ii + tcol * n + kk*ncols*n] = ci[ii].Na[jj]; tcol ++; } // get the loci information
				}
				if (kk==1){
					int tcol = 2; //define the column place holder
					for (int jj=0; jj<n_disp_loci; jj++) { out[ii + tcol * n + kk*ncols*n] = ci[ii].Db[jj]; tcol ++; } // get the loci information
					for (int jj=0; jj<n_neutral_loci; jj++) { out[ii + tcol * n + kk*ncols*n] = ci[ii].Nb[jj]; tcol ++; } // get the loci information
				}
			} 
		}
		else { // if ind dead.. set everything to -9
			for (int jj=0; jj<n*ncols*2; jj++) { out[jj] = -9; }
		}	
	}
	UNPROTECT(1);
	return ans;
} */
SEXP output_individuals(ind *ci, int n) {
        double *out; SEXP ans;
		int ncol = 2 + n_disp_loci * 2 + n_neutral_loci * 2; //define the number of columns
        PROTECT(ans = allocMatrix(REALSXP, n, ncol));
        out = REAL(ans); //pointer to output dataset
        for (int ii=0; ii<n; ii++) { //cycle through each of the individuals
                if (ci[ii].sex >=0){ //if the individual is alive, copy the necessary info
                        out[ii + 0 * n] = ci[ii].x; out[ii + 1 * n] = ci[ii].y; // set the xy values
                        int tcol = 2; //define the column place holder
                        for (int jj=0; jj<n_disp_loci; jj++) { out[ii + tcol * n] = ci[ii].Da[jj]; tcol ++; out[ii + tcol * n] = ci[ii].Db[jj]; tcol ++; } // get the loci information
                        for (int jj=0; jj<n_neutral_loci; jj++) { out[ii + tcol * n] = ci[ii].Na[jj]; tcol ++; out[ii + tcol * n] = ci[ii].Nb[jj]; tcol ++;  } // get the loci information
                } else { // if ind dead.. set everything to -9
                        for (int jj=0; jj<ncol; jj++) { out[ii + jj * n] = -9; }
                }
        }
        UNPROTECT(1);
        return ans;
}
/* function to create a matrix of density information */
SEXP output_density(int **dens, int spX, int spY) {
	double *out; SEXP ans;
	PROTECT(ans = allocMatrix(REALSXP, spY, spX));
	out = REAL(ans); //pointer to output dataset
	for (int xx=0; xx<spX; xx++) {
		for (int yy=0; yy<spY; yy++) {
			out[yy + xx * spY] = dens[xx][yy];
		}
	}
	UNPROTECT(1);
	return ans;
}


/* define the mother function
	R_n is the starting number of individuals
	R_spX & R_spY are the dimensions of the landscape
	R_d_cost is the dispersal cost */
SEXP mother(SEXP R_n, SEXP R_spX, SEXP R_spY, SEXP R_lambda, SEXP R_K, SEXP R_disp_cost, SEXP R_reports_rate, SEXP R_num_gens, SEXP rho, SEXP R_m_rate) {
	GetRNGstate();
	//define the info brought in from R
	int spX = INTEGER_VALUE(R_spX); int spY = INTEGER_VALUE(R_spY); //size of the habitat matrix
	int n = INTEGER_VALUE(R_n); //define the current population size
	int K = INTEGER_VALUE(R_K); //define the cell size carrying capacity
	double lambda = NUMERIC_VALUE(R_lambda); // define the lambda value
	double disp_cost = NUMERIC_VALUE(R_disp_cost); //define the dispersal cost
	int report_rate = INTEGER_VALUE(R_reports_rate); // get the reporting rate
	int num_gens = INTEGER_VALUE(R_num_gens); //assign the number of new generations
	double m_rate = NUMERIC_VALUE(R_m_rate); //assign mutation rate
	
	//define some other variables
	int ii, jj, xx, yy, zz; //some iteration variables
	nb_list *neighbors = Calloc(9,nb_list); initialize_neighborlist(neighbors,spX,spY); //populate and initialize the list of neighbors
	int **nbIndex = Calloc(spX,int*); for (xx=0; xx<spX; xx++) { nbIndex[xx] = Calloc(spY,int); } //allocate the memory for the nb Index values
	assign_nb_values(nbIndex,spX,spY); //get the neighbor index values
		
	//initialize the starting set of individuals
	ind * cur_ind = Calloc(n+1,ind); //define the object for the current individuals
	init_inds(cur_ind,spX,spY,n); //initialize the current indidviduals
	
	//define the initial habitat grid & populate the data
	int ** density; int ** nmales; int *** mIndex; //set the grid data holders
	density = Calloc(spX, int*); for (xx=0; xx<spX; xx++) { density[xx] = Calloc(spY,int); } //allocate the memory
	nmales = Calloc(spX, int*); for (xx=0; xx<spX; xx++) { nmales[xx] = Calloc(spY,int); } //allocate the memory
	reset_grid(density, spX, spY); reset_grid(nmales, spX, spY);// initialize the values to 0
	count_individuals(cur_ind, density, nmales, n); //populate with counts of individuals
	mIndex = Calloc(spX, int**); for (xx=0; xx<spX; xx++) { mIndex[xx] = Calloc(spY,int*);//allocate the memory
		for (yy=0; yy<spY; yy++) { mIndex[xx][yy] = Calloc(nmales[xx][yy],int); //allocate 3rd dim of this matrix
			for (zz=0; zz<nmales[xx][yy]; zz++) { mIndex[xx][yy][zz] = -99; } //set everything to -99 
		}
	}
	index_males(cur_ind, mIndex, nmales, n); // get the index values for the males in that cell
	
	//setup the "burn in" and then expand the grid and continue...
	for (ii=0; ii<500 + num_gens; ii++){ //cycle through the generations
		int new_n = 0; int *n_offspring = Calloc(n,int); //define some objects to hold temporary info
		for(jj=0; jj<n; jj++) { //define the number of offspring for the females of the current individuals
			if (cur_ind[jj].sex==1 && nmales[cur_ind[jj].x][cur_ind[jj].y]>0) { //ensure there are males to reproduce with
				double E = hass_comm(lambda, K, density[cur_ind[jj].x][cur_ind[jj].y]);
				n_offspring[jj] = Rf_rpois(E);
			} else { n_offspring[jj] = 0; } //define the number of offspring for females
			new_n += n_offspring[jj]; // keep a running tally of the new number of individuals
		}
		ind *new_ind = Calloc(new_n, ind); //define the memory for the new individuals
		reproduce(cur_ind, new_ind, n_offspring, n, nmales, mIndex, m_rate); //define the next generation individuals now that we know how many offspring each female will have
		disperse(new_ind, new_n, neighbors, nbIndex, disp_cost); // disperse individuals as necessary
	
		//convert new to current
		n = new_n; cur_ind = Realloc(cur_ind, n, ind); //allocate the memory
		copy_inds(cur_ind, new_ind, n); // copy the individuals
		
		//reextract grid information
		reset_grid(density, spX, spY); reset_grid(nmales, spX, spY);// initialize the values to 0
		count_individuals(cur_ind, density, nmales, n); //populate with counts of individuals
		for (xx=0; xx<spX; xx++) {
			for (yy=0; yy<spY; yy++) { 
				if (nmales[xx][yy]==0) { mIndex[xx][yy] = Realloc(mIndex[xx][yy],1,int); 
				} else { mIndex[xx][yy] = Realloc(mIndex[xx][yy],nmales[xx][yy],int); }//allocate 3rd dim of this matrix
				for (zz=0; zz<nmales[xx][yy]; zz++) { mIndex[xx][yy][zz] = -99; } //set everything to -99 
			}
		}
		index_males(cur_ind, mIndex, nmales, n); // get the index values for the males in that cell
		Free(n_offspring); Free(new_ind); //free up the memory
		
		//if the burn in time is complete, expand our grids to allow further spread
		if (ii == 500) { printf("burn in done and current population size is %i\n", n);
			for (xx=0; xx<spX; xx++) { Free(density[xx]); Free(nmales[xx]); Free(nbIndex[xx]); }
			for (xx=0; xx<spX; xx++) { for (yy=0; yy<spY; yy++) { Free(mIndex[xx][yy]); } Free(mIndex[xx]); }
			spX+=num_gens; //increase spX size
			nbIndex = Realloc(nbIndex,spX,int*); for (xx=0; xx<spX; xx++) { nbIndex[xx] = Calloc(spY,int); } //allocate the memory for the nb Index values
			assign_nb_values(nbIndex,spX,spY); //get the neighbor index values
			density = Realloc(density, spX, int*); for (xx=0; xx<spX; xx++) { density[xx] = Calloc(spY,int); } //allocate the memory
			nmales = Realloc(nmales, spX, int*); for (xx=0; xx<spX; xx++) { nmales[xx] = Calloc(spY,int); } //allocate the memory
			reset_grid(density, spX, spY); reset_grid(nmales, spX, spY);// initialize the values to 0
			count_individuals(cur_ind, density, nmales, n); //populate with counts of individuals
			mIndex = Realloc(mIndex, spX, int**); for (xx=0; xx<spX; xx++) { mIndex[xx] = Calloc(spY,int*);//allocate the memory
				for (yy=0; yy<spY; yy++) { mIndex[xx][yy] = Calloc(nmales[xx][yy],int); //allocate 3rd dim of this matrix
					for (zz=0; zz<nmales[xx][yy]; zz++) { mIndex[xx][yy][zz] = -99; } //set everything to -99 
				}
			}
			index_males(cur_ind, mIndex, nmales, n); // get the index values for the males in that cell
		}
		// store key information out into R at specific times
		if (ii>=500) { //if past the burn in time
			if ( (ii%500) % report_rate == 0 || ii == 500 + num_gens - 1) { // if it is the last gen or at the expected reporting rate
				char tname[10]; sprintf(tname, "%s%d", "gen_", ii-500); //setup the dynamic name for R
				defineVar(install(tname), output_individuals(cur_ind,n), rho); //write a matrix of individuals out to R
				sprintf(tname, "%s%d", "dens_", ii-500); //setup the dynamic name for R
				defineVar(install(tname), output_density(density,spX,spY), rho); //write a matrix of individuals out to R
			}
		}
	}
	printf("Final population size is %i\n", n);
	/*
	//perform some checks
	check_nb_values(nbIndex,spX,spY); Rprintf("\n"); //check the neighbor zone values
	check_inds (cur_ind, spX, spY, 10); Rprintf("\n"); //check the first 10 individuals
	check_grid(density, spX, spY); Rprintf("\n"); //check the values
	check_grid(nmales, spX, spY); Rprintf("\n"); //check the values
	check_mIndex (cur_ind, mIndex, nmales,1,1); Rprintf("\n"); //check the indexing
	*/
	
	//return some useless info
	SEXP res; PROTECT(res = allocVector(INTSXP,1));
	INTEGER(res)[0] = -9999;

	PutRNGstate();
	//release and free up memory
	Free(neighbors); Free(cur_ind);
	for (xx=0; xx<spX; xx++) { Free(density[xx]); Free(nmales[xx]); Free(nbIndex[xx]); }; Free(density); Free(nmales); Free(nbIndex);
	for (xx=0; xx<spX; xx++) { for (yy=0; yy<spY; yy++) { Free(mIndex[xx][yy]); } Free(mIndex[xx]); } Free(mIndex);
	
	
	UNPROTECT(1);
    return res;
}

